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2.5D  geometry  and  is  expounded  in  detail  in  a  recent  Center  for  Wave 
Phenomena  Research  Report  (Bleistein,  1984a). 

Beylkin's  paper  uses  powerful  mathematical  tools,  such  as  the  notions 
of  pseudo-differential  operators,  generalized  Radon  transforms,  and 
generalized  back  projections.  Moreover,  Beylkin  frames  his  work  in  an  N- 
dimensional  space. 

Here,  we  dispense  with  this  mathematical  machinery  and  for  convenience 
confine  ourselves  to  the  3D  case  (and  its  2.5D  specialization).  We  are  able 
to  expound  Beylkin's  results  by  an  approach  similar  to  that  presented 
earlier  in  (Cohen  and  Hagin,  1985).  However,  we  do  not  attempt  to 
rigorously  prove  our  results,  but  instead  content  ourselves  with  an 
intuitive  derivation  and  the  citing  of  Beylkin's  theorem.  We  remark  in 
passing  that  we  are  contemplating  the  possibility  of  developing  a  classical 
proof  of  Beylkin's  result. 

We  have  established  significant  simplifications  of  Beylkin’s  result 
for  the  following  cases  of  propagation  governed  by  the  acoustic  wave 
eq  uat ion : 

(1)  The  zero-offset  case  with  a  general  c(x,y,z)  reference  velocity.  This 
work  is  described  in  the  present  report.  See  also  (Cohen  and  Hagin, 
1985)  for  the  special  case  of  a  c(z)  reference  velocity.  Here,  in  all 
but  the  case  of  a  constant  reference,  one  has  to  exclude  "turned  rays" 
in  order  to  guarantee  that  h  (or  its  equivalent  in  earlier  approaches) 
does  not  vanish.  It  is  easy  to  build  this  restriction  into  the  zero- 
offset  algorithm  (Sumner,  1985),  but  only  at  the  cost  of  excluding  the 
imaging  of  reflectors  "from  below".  At  the  present  time,  the  Center 
for  Wave  Phenomena  is  investigating  the  use  of  turned  rays  in  the  c(z) 


(2)  The  common  source  offset  gather  and  the  common  receiver  offset  gather 
with  general  c(x,y, z)  reference  velocity.  This  work  is  also  described 
below . 

(3)  The  case  of  common  offset  with  a  constant  reference  velocity.  Here,  an 
explicit  inversion  formula  has  been  obtained  (Snllivan  and  Cohen, 
1985).  Ve  remark  that  the  derivation  of  this  result  does  not  proceed 
from  the  Born  approach  described  below,  but  instead  starts  from  the 
Kirchhoff  (high  frequency)  representation  of  the  scattered  field 
(Bleistein,  1984b).  Work  on  coding  and  testing  this  algorithm  is  now 
under  way. 


IIi_  The  General  Imweraion  Formula 

We  consider  a  completely  arbitrary  source-receiver  configuration 
parametrized  by  two  surface  coordinates,  and  Vectors,  x  and  x' 

denote  arbitrary  subsurface  field  points.  Vector  jr  denotes  a  generic 
source,  while  xfi  denotes  a  generic  receiver.  Any  relation  between  the 
sources  and  receivers  is  defined  by  the  dependencies  of  x  and  x  on  the 
parameters  and  ?a: 


x  =  <wV 


£  =  <«*.«,> 


z 

-s 


X  9  -  (x'a  X',  X') 

-  i  a  s 

=xs(£).  5r=xr(i>. 


(i) 


For  example,  if 

;r<£>  -  st(£)  . 


(2) 


then  we  have  the  case  of  zero  offset.  If 

x  ( l )  =  constant, 

“  8  * 


(3) 


then  we  have  the  common  source  configuration 


If 


xt(£)  =  xs(£)  +  2d,  d  =  constant. 


(4) 


then  we  have  the  common  offset  case,  etc. 

Ve  note  that  the  above  foronlation  of  the  sonrce-receiver  configuration 
is  somewhat  more  general  and  inclusive  than  the  formulation  in  Beylkin’s 
paper.  Our  formulation  admits  the  possibility  of  a  curved  observation 
surface,  while  he  makes  the  usual  assumption  of  observations  on  the  flat 
observation  plane  z  =  x,  =  0.  Ve  have  included  this  mild  generalization 
because  it  may  aid  in  the  development  of  pre-statics  inversion  algorithms 
(May  and  Covey,  1980).  Also,  Beylkin  splits  off  his  derivation  of  the 
common  source  result  into  a  separate  case  from  the  case  when  the  receiver 
locations  depend  on  the  source  locations. 

As  expounded  in  several  previous  papers  and  reports,  we  make  the 
following  assumptions: 

(1)  the  velocity.  v(x,y, z),  is  well  approximated  by  a  known  reference 
velocity,  c(x,y,  z),  so  that 


— -  =  — j -  ( 1  +a(  x) ) ,  (5) 

v  (x)  c  (x) 

where  a(x,y, z)  is  a  perturbation  correction, 

(2)  the  seismic  fields  are  governed  with  sufficient  accuracy  by  the  3-D 
acoustic  wave  equation,  and 


(3)  the  seismic  source  can  be  reduced  to  an  ideal  3-D  point  source,  so 


These  assumptions  and  the  application  of  Green's  theorem  lead  to  the  linear 
integral  equation. 


tc[  G  ( w,  x ' ,  x  )  G  ( to,  x ' ,  x  ) 

D(w,xr,xs)  =  w  jJJ  d  x'  - * - —  Q(x') 


c*(x') 


(7) 


Here,  by  the  assumption  that  a  in  (5)  is  small,  we  have  been  able  to 

approximate  the  field  emanating  from  the  sources  by  one  governed  by  a  wave 
equation  with  a  velocity  function  equal  to  the  reference  velocity,  c, 

instead  the  full  velocity,  v.  Since  our  sources  are  point  sources,  the 

functions  denoted  by  G  in  equation  (7)  are  the  impulse  response  or  Green's 

functions  for  the  wave  equation  with  velocity  c(x,y, z).  Equation  (7)  is  an 
integral  equation  for  the  unknown  velocity  perturbation,  a(x,y, z),  with  data 
being  the  observations,  D,  at  the  receivers,  xr>  due  to  the  sources,  x$. 

We  now  exploit  the  fact  that  geophysical  data  resides  in  the  high 
frequency  regime  (Bleistein,  1984b),  so  that  we  may  replace  the  Green’s 
functions  in  (7),  by  their  WKBJ  approximations, 


G(fa>,  x,  x  ) 

-  -o 


A(x. xo) 


iwr(x,x  ) 
e  -  “o 


(8) 


where  xQ  is  either  xs  or  xr  and  the  travel-time  phase  satisfies  the  eikonal 
equation. 


Vx-Vv  =  — -  . 

c  (x) 

while  the  amplitude  satisfies  the  transport  equation. 


(9) 


2Vt-VA  +  AV*t  =  0 


(10) 


Thus,  we  may  rewrite  (7)  as 


D(<o,  x  ,  x  ) 
—  r  s 


q(  X  ') 


(11) 


'  .  *  -  c  '-'f-'T-  ■  -  »-»I'  T  "•  ».» 


.*  .  a(x',£)  _  iu>l(  x '»  $) 
o  x  “ — —  e  " 

c  (x') 

where  |(x,£)  and  a(x,£)  are  given  by 

l(x.£)  =  r(x,x$)  +  r(x.xr)  , 
a(x,£)  =  A(x,  xg)  A(x,  xf) 


(12) 


From  the  basic  migration  principle  of  "backward  propagation",  it  is  not 
hard  to  gness  that  the  inversion  operator  for  (11)  will  have  the  negative  of 
the  phase  in  (11).  The  correct  amplitude  is  not  as  easy  to  gness  and  in  onr 
inversion  equation, 

a(x)  -  jj  d’$  jiTIT  b(x.£)  |  du>  FU)  e_1U,(?' ^  D(w,  xf,  xs)  ,  (13) 


we  merely  denote  it  by  b(x,£),  and  deduce  it  below.  In  equation  (13),  (|g 
denotes  the  differential  area  element  on  the  observation  surface  (unity 
for  a  plane)  and  F(u)  denotes  a  known  high  band  pass  filter  which  is 
included  to  honor  the  fact  that  the  data  is  confined  to  the  high  frequency 
regime.  Note  that  while  a  priori  the  unknown  amplitude  conld  depend  on  w, 
we  suppress  this  potential  dependence  because  our  prior  inversion  results 
suggest  that  b(x, £)  is  independent  of  u  to  first  order.  The  results  to 
follow  also  confirms  this  assumption. 


Inserting  (11)  into  (13)  gives  an  equation  which  maps  a(x’)  to  a(x)  . 
Thus,  since  we  have  an  integral  over  all  x'  space,  the  kernel  must  be  the 
three-dimensional  Dirac  delta  function.  A  slight  rewriting  of  this  fact 


f|  <3*5  Jg  a(x'.£)b(x.£)  • 

f  .  *17/  »  iwfi(x', £) -|(x, £) ]  *,  ...  M.v 

I  din  in  F(w)e  -  c  ( x ) 6 ( x-x  )  .  (14) 

Oor  task  is  to  solve  (14)  fot  the  unknown  amplitude,  b.  Then  equation  (13) 
will  give  the  inversion  algorithm  for  a. 

Since  we  have  used  high  frequency  approximations,  we  cannot  recover  low 
frequency  phenomena  such  as  trend  information.  We  can  expect  to  recover 
information  about  the  discontinuities  (i.e.  rapid  changes)  in  the 
substructure . 

The  essence  of  Beylkin's  result  is  that  if  we  confine  attention  to  the 
discontinuity  structure  of  a,  then  in  (14),  we  need  keep  only  two  terms  of 
the  Taylor  series  of  |(x',£)  about  x ’  =  xr  and  similarly  we  need  keep  only 
one  term  of  the  expansion  of  »(*',$).  In  order  to  establish  these  facts 
rigorously  Beylkin  uses  powerful  results  in  the  theory  of  generalized  Radon 
transforms.  However,  noting  that  the  result  on  the  left  must  involve  a 
Dirac  delta  function  acting  at  x'  =  x»  these  approximations  are  intuitively 
reasonable.  Thus  in  (14),  we  use  the  approximations 

!(x',£)  *  Kx,£)  +  V|(x,  ^)  •  (x'-x) 

•  (x  ',£)  «  a (x, £)  (15) 

to  obtain 


||  d*£  ^g  a(x,£)b(x,£)  * 

|  dw  o)aF(u»)eih>Vf b  c* (x) M x-x ' )  .  (16) 

Following  Beylkin,  we  make  the  change  of  variables, 

k  =  <dV|(x,£)  (17) 

from  <«>.  tj,  to  kt,  ka,  k,  (while  viewing  x  as  a  'parameter").  Since  the 


to  dependence  on  the  Jacobian  of  this  transformation  is  jnst  a  power,  Beylkin 
defines  the  scaled  Jacobian,  h,  by 


31 

31 

a* 

3x 

3x 

dx 

i 

a 

s 

„a_ 

,  a 

z 

3  1 

3V$a 

5_| - 

0*a3$a 

3*1 

3*1 

3*1 

av>sa 

3*a35* 

We  now  rewrite  (16)  as 

/II  d>k  Jh TxT'£)~|  ~  }  *  c*(x)6(x-x' 


From  (19)  and  the  classical  Fonrier  transform  completeness 
is  clear  how  to  choose  the  inversion  amplitude  b: 


b(x,£) 


1  c*(x)lb(^>l 

8n*  ^g  a(x,£) 


Inserting  this  result  into  (13)  gives  the  inversion  result: 

•<«>  -  //  A 

where  |  and  a  are  defined  above  in  (12). 

It  is  intriguing  to  note  that  if  we  define  the  kernel. 


K(  to,  x,  i  ,  x  )  = 
-  -s  -r 


<u  G(w,  x,xs)G(w,  x,  xf) 


a . 

c  (x) 


(18) 


(19) 


relation,  it 


(20) 


(21) 


then  we  may  state  our  integral  equation  and  its  inversion  as 


III^_  The  Zeto-Offaet  Conf iguration 

In  the  case  of  zero-offset,  we  can  find  a  representation  of  the 
Jacobian  h  in  terms  of  the  rays  used  to  construct  the  Green's  function. 
Using  the  notation  x0>  as  in  (2),  for  the  zero-offset  source-receiver  point 
we  have 


2Vt(x,x  ) 

-  “0 

E< 

s'?o) 

h(x,£)  = 

2  Vt  ( x,  x  ) 

2<U  ■  0 

=  8 

a 

3^ 

p(x,  X  ) 

K  -  0 

2  Vt  ( x,  x  ) 

- 0 

a 

HI 

HI 

CU 

where  in  the  second  equality,  we  have  introduced 

p(x,  x  )  =  Vt(x,  x  )  . 

—  -  -o  -  -o 

and  where  the  £  dependence  in  p  comes  from  x0  =  xQ(£). 
From  the  eikonal  equation  (9),  we  have 


(25) 


(26) 


E 


•  p  = 


c  ( x) 


(27) 


Since  the  speed,  c,  is  independent  of  J;,  we  also  derive  from  (27) 


E  ‘ 


i  =  1,2. 


(28) 


We  now  multiply  the  third  column  of  h  by  p,  and  compensate  with  the 
reciprocal  of  p,  outside  the  determinant/  then  on  multiplying  the  first  two 
columns  by  respectively  p2  and  p2*  and  adding  these  to  the  third  column,  we 


ob  t  a  i  n 


Mx,fc>  =  -2- 


pl 

pa 

3Pi 

3pa 

05a 

05a 

api 

3pa 

3*a 

35. 

which  can  be  written  as 


g  a(Pi'P*) 

*  a  a (?  ,5  ) 

c  (5)p,U»s0>  ?1 


Here  we  have  introduced  the  Jacobian  notation 


dtPi-  Pa) 


3pi 

3Pj 

3?a 

05a 

3pi 

3pa 

3*a 

3?a 

The  evaluation  of  this  determinant  is  a  key,  and  often  difficult,  step.  In 
appendix  A,  we  show  that 

3(px'p*)  ,  a 

=  '16n  A  (5'5o)p»{5'50)  8  -E0  (32) 

where  g0  is  the  initial  direction  of  the  ray  from  xQ  to  x  and  n  is  the 

downward  normal  to  the  data  surface  at  x  .  Thus  the  absolute  value  of  h  is 

-0 


given  by 


8A*<*,x> 

|h(*.ft)  I  =  — -  '  Jg  8  *E0 

c  (x)  1 


and  the  zero-offset  inversion  for  general  reference  velocity  can  be 
expressed  as 


-  11 


~  ff  d*5  J8  S  *  Eo  f  dw 


-2iwr(x,  x  ) 


D  ( <D,  x  )  • 
-o 


This  is  >  substantial  simplification  of  the  general  result  (21),  since  it 
obviates  the  necessity  for  computing  the  amplitude  when  tracing  the  rays 
between  surface  and  field  points.  Note  from  appendix  A  that  computation  of 
the  amplitude  would,  in  turn,  require  computation  of  the  ray  Jacobian.  The 
net  result  is  that  we  save  six  equations  along  rays  and  have  only  to  solve 
the  remaining  seven  equations  for  x,  p>  and  z. 

The  special  case  of  this  result  in  which  the  reference  speed  is  a 
function  of  depth  only  was  derived  in  (Cohen  and  Hagin  1985).  That  result 
agrees  with  the  present  one  aside  from  a  slight  change  in  notation.  Of 
course  the  present  result  also  agrees  with  the  constant  reference  speed 
result  (Bleistein,  1985). 

There  are  two  further  issues  to  be  addressed  concerning  the 
implementation  of  (34).  These  are  the  specialization  to  the  case  of  a 
linear  array  (the  2.5  dimensional  case  mentioned  in  the  introduction)  and 
the  fact  that  discerning  discontinuities  is  made  easier  if  one  replaces  a  by 
the  reflectivity  function,  p,  (related  to  the  normal  derivative  of  a) 
(Bleistein,  1984b). 

In  the  usual  case  of  a  data  set  collected  on  a  linear  array  instead  of 
a  full  20  array,  we  cannot  obtain  a  full  3D  subsurface  inversion.  A  model 
consistent  with  this  restricted  set  of  observations  is  one  in  which  both  the 
observation  surface  and  the  subsurface  are  assumed  to  be  "cylindrical’'  with 
respect  to  x, : 

c  =  c(x  ,x  )  ,  a  =  a(x  ,x  )  (35) 


k-  -  ■-*1 


.1.,.  »  ■  v  ■  i 


and  that  the  observation  surface  had  the  special  parametrization: 


‘x=  Xx<U' 


‘»=  *»' 


x  =x  (5  ) 

1  »  ’1 


(36) 


In  the  2.5  case,  it  is  possible  to  elioinate  the  integral  over  £a  by  the 
method  of  stationary  phase.  Carrying  ont  this  step  and  for  simplicity 
further  specializing  to  the  case  of  linear  observations. 


we  obtain  the  2.5D  result. 


a(xi,x|) 


-2iut(x  ,  x  ,  £  ,  0) 

D(o).^) 


(3  8) 


Here  af  denotes  the  final  value  of  the  ray  parameter  o  for  the  ray  from  x0 
to  x  (or  equivalently  from  x  to  x0)  and  q,  is  the  valne  of  pg  on  the 
observation  array  (see  (40)  below).  In  the  reference  just  cited,  it  is 
shown  that  we  may  use  2D  calculations  to  compute  the  rays: 


dxi 

JT  ~  P1  '  V0)  =  *1 


x^O) 


0 


dPj 

do  7  t>x 


c  (*a»*.) 


PJO)  -  q4 


t(0)  =  0 


(39) 


To  apply  (38).  we  fix  the  desired  field  point,  (xx, x, )  and  a  suitable 
aperture  of  observation  points  ($1(0)  on  the  observation  array.  Then  we 
trace  the  ray  from  ($i»0)  to  (xx, xg)  and  compute  t.  Of  and  the  ray  parameter 


Then  we  can  explicitly  evaluate 


Since  the  last  two  rows  of  each  of  these  determinants  is  the  same  as  in  the 


zero-offset  case  (with  ;r  taking  the  place  of  xQ),  we  can  proceed  as  at  the 
beginning  of  the  last  section  to  obtain 


=  P,(?**r 


J7TT  psV  V+  V 


®(P1»  Pa} 


and 


h  = 
2 


atPj.  p8) 


P4 (x, xf) c* (x)  a(W 


(48) 


(49) 


Dsing  the  resnlt  of  the  appendix,  we  have 


P  P/  p'  ♦  p’  ♦ 

»  1  »  c  (x)J 


(50) 


and  so,  nsing  the  general  result  (21)  we  obtain 


[?)  =  iff d*5  K  Mirrr  1  +  “ r1  pS|(prps >  8r  * 

$  *  3 


(  -iwfrfx.x  )  +  t(x,x  )] 

•  J  do)  F(w)  e  r 


D(w,  ;r)  • 


In  (51)  ve  have  restated  the  resolt  in  (50)  by  nse  of  the  eikonal 
equation  for  j>s> 

The  implementation  of  (51)  is  not  as  simple  as  the  zero-offset  case/ 
the  principle  difficulty  being  the  necessity  to  compote  the  two  amplitude 
factors.  In  addition,  the  questions  of  constructing  'synthetic  apertures* 
in  order  to  use  an  array  of  common  source  gathers  as  well  as  how  to 
systematically  use  nearby  conanon  source  gathers  to  reduce  noise  in  a  post 
migration  processing  remain  open. 

Obviously,  the  solution  for  the  common  receiver  configuration  can  be 
derived  at  once  from  (51)  by  interchanging  the  subscripts  s  and  r  and 


defining 


g  =  Vr(x,  x#) 


p_=  Vr ( x, x  ) 
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A1 :  Appendix  on  the  Zero-Offset  Jacobinn 

It  is  natural  to  think  of  t(x,  x0)  as  the  travel  time  along  a  ray  path 
with  initial  point  xQ  and  running  point  x.  In  this  model,  p  is  the 

associated  ray  direction  along  the  ray  from  xQ  to  x •  Furthermore,  to  solve 
the  eikonal  equation  in  this  model  by  the  method  of  characteristics  (or  "ray 
method”),  one  most  consider  the  family  of  rays  emanating  from  the  'fixed* 
point  xQ  and  consider  x  as  an  arbitrary  field  point.  Of  course,  in 

constructing  the  eikonal,  t,  and  transport  amplitude.  A,  we  can  consider  xD 

to  be  an  arbitrary  point  on  the  observation  surface,  or  otherwise  put,  we 

repeat  the  ray  method  solution  for  each  xQ. 

On  the  other  hand,  if  we  think  of  the  field  point,  x,  as  fixed,  it  is 
equally  valid  to  compute  the  travel  time  using  a  ray  with  initial  point  x 
and  ”final”  point,  xq.  In  this  model,  this  "upward”  ray  direction  is 

simply  the  negative  of  the  ray  direction  on  the  "downward"  ray  at  any  given 
point  on  the  ray,  and  we  consider  the  family  of  rays  emanating  from  the 

"fixed”  point  x. 

In  the  computation  of  the  2  by  2  Jacobian  determinant  in  (31),  we  need 
to  compute  derivatives  with  respect  to  the  second  argument  of  g.  Urns  we 
adopt  the  point  of  view  where  x  is  fixed.  However,  in  our  original 

description,  the  vector  x^  was  bound  to  the  observation  surface.  Thus  we 

introduce  a  new  running  variable,  jr,  governed  by  the  ray  eqnations. 
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We  impose  no  condition  on  the  starting  ray  direction  in  (A-l)  because,  to 


construct  the  Green's  function,  we  want  the  "conoidal”  solution  of  the 
eikonal  equation  (Bleistein,  1984b).  In  order  to  impose  the  condition  that 
each  ray  emanating  from  x  passes  through  a  point,  x0»  on  the  observation 
surface,  we  demand 

y(a{)  =  xq(£)  ( A-2) 


This  last  condition  defines,  for  each  such  ray  the  final  value,  Of>  of  o  as 
a  function  of  £.  (The  question  of  having  a  unique  ray  to  each  surface  point 
is  subsumed  in  the  issue  of  the  non-vanishing  of  Beylkin’s  b  determinant.) 

We  now  evaluate  our  2  by  2  Jacobian  by  introducing  the  ray  parameter  o 
and  then  using  the  chain  rule  as  follows: 


3  (p4  (x,  xq)  ,  p2  (x,  xq)  ) 


3(p  . p  ,  a) 
1  2 

au^.a) 


a(px  z)  *  p,  (*»  z)) 


0-0. 
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3  (  p  ,  p  ,  o) 
1  2 


9(y1.yi,yt) 


a(yj.yt.y,) 


a  =of  9($x.  5#,  o) 


a  =  a. 


The  second  of  these  3  by  3  Jacobians  can  be  evaluated  as 


a(yi'y*’yi) 


is 

0=0- 

t. 

f 

-1 

■B0  *  ‘x  *  »» 


-to  •  * 


where  we  have  introduced  the  surface  tangent  vectors 


( A-4) 


-i  ”  9^ .  -o  '  i  =  1*  2  '  (A-5) 

and  where  j>0  denotes  the  initial  direction  of  the  ray  from  x0  to  x  (or  the 
negative  of  the  final  direction  of  the  ray  from  x  to  xQ)  and  finally  8 

denotes  the  nnit  downward  normal  to  the  data  surface  at  xq.  The  remaining 

Jacobian  in  (A-3)  is  the  reciprocal  of  the  ray  Jacobian  for  rays  emanating 
from  x 

a(yi'yi'y») 

J(y)  =  atp^.p’.o)  *  ( A-6 ) 

From  the  transport  equation,  we  can  derive  the  relation, 

J(y)  A  (y, x)  =  constant  ( A — 7 ) 

(see  (Bleistein,  1984b)).  The  constant  in  ( A— 7 )  can  be  evaluated  by 

allowing  j  to  approach  the  fixed  field  point  x  and  using  the  constant  c 

result  with  c  being  c(x)<  We  find  for  this  limit. 


lim  JA* 

z-»* 


1 _ 

16jt  p^  (x,  xq) 


( A-8) 


so  that  the  ray  Jacobian  is  given  by 


J 


=  o- 


( A-9) 


Here  we  have  used  reciprocity  to  switch  the  arguments  of  A.  Combining  ( A- 
3),  (A-4),  and  (A-9)  we  find  that  the  required  2  by  2  Jacobian  is  given  by 


9( p  ,  p  ) 
»  a 


( A-10) 


as  asserted  in  the  text  (32). 
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